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ABSTRACT 

The long-term evolution of massive black hole binaries at the centers of galaxies is studied in a variety 
of physical regimes, with the aim of resolving the "final parsec problem," i.e. how black hole binaries 
manage to shrink to separations at which emission of gravity waves becomes efficient. A binary ejects 
stars by the gravitational slingshot and carves out a loss cone in the host galaxy. Continued decay 
of the binary requires a refilling of the loss cone. We show that the standard treatment of loss cone 
refilling, derived for collisionally relaxed systems like globular clusters, can substantially underestimate 
the refilling rates in galactic nuclei. We derive expressions for non-equilibrium loss-cone dynamics and 
calculate time scales for the decay of massive black hole binaries following galaxy mergers, obtaining 
significantly higher decay rates than heretofore. Even in the absence of two-body relaxation, decay 
of binaries can persist due to repeated ejection of stars returning to the nucleus on eccentric orbits. 
We show that this recycling of stars leads to a gradual, approximately logarithmic dependence of the 
binary binding energy on time. We derive an expression for the loss cone refilling induced by the 
Brownian motion of a black hole binary. We also show that numerical iV-body experiments are not 
well suited to probe these mechanisms over long times due to spurious relaxation. 
Subject headings: black hole physics — galaxies: nuclei - stellar dynamics 



1. INTRODUCTION 

Larger galaxies grow through the agglomeration of 
smaller galaxies and protogalactic fragments. If more 
than one of the fragments contains a massive black 
hole (MBH), the MB Hs will form a bound system in 
the merger product llBeeelman. Blandford fc Reeslll98(i 
lRoollT98l iValtaoia. Valtonen fc Bvrdlll989f) . This sce- 
nario has received considerable attention because the ul- 
timate coalescence of such a pair would generate an ob- 
servab le outburst of gravity waves llThorne fc Bra'g inskii 
HUH). iBegelman. Blandford fc Reesl (jl980f ) showed that 
the evolution of a MBH binary can be divided into three 
distinct phases: 1. As the galaxies merge, the MBHs 
sink to ward the center of the new galaxy via dynamical 
friction l)Chandrasekharll943[) where they form a binary. 
2. The binary con tinues to decay via gravitational sli ng- 
shot interactions iSaslaw. Valtonen fc Aarsethl 1197-3) in 
which stars on orbits intersecting the binary's are ejected 
at velocities comparable to the binary's orbital velocity, 
while the binary's binding energy increases. 3. If the bi- 
nary's separation decreases to the point where the emis- 
sion of gravity waves becomes efficient at carrying away 
the last remaining angular momentum, the MBHs coa- 
lesce rapidly. In this paper we explore Phase 2 and its 
transition to Phase 3; this transition is understood to be 
the bottleneck of a MBH binary's path to coalescen ce. 

The picture of lBegelman. B landfor d fc Reesl l)1980j) has 
remained essentially valid in the face of several impor- 
tant developments in intervening years. Early theoretical 
studies envisioned the centers of galaxies as constant- 
density cores, with small, embedded stellar cusps that 
develop around centra l MBHs via collisional relaxation 
iBahcall fc Wolflll97fl ll^sellT97cl IPeebleslll972t lYounel 
I1977T) . Observations at the time lacked the resolution to 
verify this picture. In the past decade, space-based obser- 



vations have revealed that, in the majority of early-type 
galaxies, the central luminosity density increases contin- 
ually inward as an approximate power law p(r) ~ r 1 
down to regi ons where the MBH dominates the gravita- 
tional force llFerrarese et alJll994t iCebhardt et aTlll99fit 
iRest et al.l200lh . The tvnical luminosity profile of a faint 
elliptical galaxy is a nearly-featureless power law with 
1.5 < 7 < 2.5 at the resolution limit. The times required 
to build cusps via two-body relaxation are longer than 
the a ge of the Universe in most nuclei (e.g. iFaber et alJ 
Il997[) and so relaxation is not a viable explanation for 
how the cusps formed. In this work we assume that 
galaxies inherited their steep central density profiles from 
an early epoch in which the MBH grew. Cold dark 
matter (CDM) halos produced by hierarchical struc- 
ture formation that presently host galaxi es and MBHs 
were cuspy at the time o f virialization (e.g. lGhigna et all 
l2000t IPower et al]l2003|) . On smaller scales, the stellar 
cusps may have resulted from gaseous dissipation of star- 
forming medium inside pre-existing CDM cusps, from 
adiabatic response of the stars to growth of the MBH, 
or both. The same processes that led to cusp forma- 
tion were likely to have fueled the early growth of MBHs 
and, thus, set the MBHs and galaxies on a course of co- 
evolution. We therefore assume that the loci of MBHs 
have coincided with the centers of cosmic overdensities 
(dark matter and stellar cusps) since the time MBHs first 
ascended into the massive range, M, ~ 1O 6 M . 

The cusps that accompany the MBHs endow 
them with an enlarged ef fective dynamical mass. 
Milos avTievic fc MerrittJ 1)200 If) found that the dynamical 
friction time scale (Phase 1) is a function of this effective 
mass and much shorter than the orbital decay time im- 
plied by the masses of the MBHs alone. This is of critical 
importance for the capacity of MBHs originating in dif- 
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ferent, merging galaxies to capture each other and form 
a hard binary in the resulting galaxy. A MBH binary is 
customarily called "hard" if its binding energy per unit 
mass exceeds the specific kinetic energy of ambient stars, 
Gp/Aa > a 2 , where jU = M 1 M 2 /(M 1 +M 2 ) is the reduced 
mass, a is the binary's semi- major axis and a is the 
ID st ellar velocity dispersion. IMilosavlievic fe MerrittJ 
( 2001) also found that in mergers of equal-mass p(r) ~ 
r~ 2 galaxies, the true separations at which the binaries 
settle at the beginning of Phase 2 are ~ 10 times smaller 
than implied by arguments based on the equipartition of 
energy. 

Masses of MBH's correlate exceptionally well with 
the central l ine-of-sight velocity dispersions of their 
host galaxies l|Ferrarese fc Merrittll2000T: iGebhardt et alJ 
l2000jk It is likely that the masses are also related 
to the depth of the potential wells associated with 
the dark matter halos in w hich the host galaxies are 
embedded ijFerrare sc 2002a), which justifies attempts to 
characterize the evolution of the MBH population using 
the statistics of hierarchical fo rmation of structure in 
dark matt er cosmologies ([ Kauffm ann fc HaehneltJ 
20001 iMenou. Haiman fc Narayanan! 120011: 

Volonteri, Ha ardt fc Madaul l2002[) . As a corollary, 
one can estimate the frequency of MBH binary coales- 
cence at redshifts that are beyond the reach of existing 
telescopes but within reach of a future space-based 
gravity wave detector such as LISA. 1 The coalescence of 
two 10 5 — 10 6 Af Q MBHs occurring practically anywhere 
in the observable Universe would generate a significant 
signal for LISA. 

Early calculations (e.g. iHaehnelti Il994(l identified the 
MBH coalescence event rate with the galaxy merger 
rate. It was subsequently noted that stellar-dynamical 
processes are inefficient at facilitating the decay of 
MBH binaries when their separations shrink below 
about one parsec and that the ultimate coalescence 
should not be taken f or granted llGquld fc Rid 120001: 
IMilosavlievic fc Merrittil2lifmiValtonedll996atlYiJl20f^ 
Circumstantial evidence, however, suggests that long- 
lived MBH binaries are in fact rare. The coincidence of 
mean black hole mass densities derived from kinematical 
studies of local galaxies and distant quasars l)Ferraresel 
(2002b)) suggests that MBHs are rarely ejected from 
galaxy centers by the strong interactions that would ac- 
company the infall of a third MBH into a nucleus contain- 
ing an uncoalesced binary. Furthermore there is evidence 
for MBH coalescences in the morphology of radio galax- 
ies, a t a rate approximately equal to the galaxy merger 
rate ijMerritt fc Ekersll20fT2]) . 

There are many processes that could drive the MBH 
binaries that form in galaxy mergers to coalescence. In 
this paper we present some new ideas regarding one such 
process, the gravitational slingshot ejection of stars. 

Ejection of stars from the galaxy nucleus will result 
in a depletion of stars in the "loss cone." We define 
the loss cone as the phase space domain in which or- 
bits approach the binary closely enough to be influenced 
by the near field of the individual binary components. 
The orbits in the loss cone are subject to gravitational 
capture and slingshot ejection. The loss cone has thus 
traditionally been viewed as a sink in the phase space, 

1 Laser Interferometer Space Antenna, |http:/ /lisa.jpl.nasa.gov| 



surrounded by a distribution of stars in equilibrium with 
respect to two-body gravitational encounters. The latter 
assumption allows a steady-state solution for the phase 
spac e density near the loss c one boundary to be derived 
fe.g. iCohn fc Kulsrudlll978lk this steady state reflects a 
balance between the supply of stars to the loss cone by 
gravitational scattering and the loss due to capture or 
ejection. In this standard picture, stars that scatter into 
the loss cone are ejected from the system with net energy 
loss, resulting in a more tightly bound binary. After all 
of the stars originally inside the loss cone have been re- 
moved, continued decay of the binary hinges on the flux 
of stars scattering into the loss cone. 

In this paper we study the physics of the "worst case 
scenario," involving nearly spherical galaxies, in which 
the MBH binaries are least likely to coalesce in their life- 
time. In §|21 we review the standard solution for the flux. 
In § we question the applicability of some of the as- 
sumptions that were made to derive that solution. In § 0] 
and §0we introduce some novel aspects of the dynamics 
of stars in the loss cone and outline an iterative proce- 
dure for calculating the long-term evolution of a massive 
MBH binary in a spherical galaxy. The details empha- 
sized here tend to shorten the timescale for the decay of 
MBH binaries, sometimes significantly. In § H3 we eval- 
uate the effectiveness of A r -body simulations to model 
the long-term evolution of MBH binaries. In §0we ad- 
dress the implications of the Brownian motion of MBH 
binaries for the long-term evolution. In § |HJwe review 
the likelihood for the coalescence of MBH binaries. Ex- 
tensions to axisymmctric and triaxial galaxies, as well as 
applications to specific galaxy models, are deferred for 
future work. 

2. REVIEW OF STANDARD LOSS CONE THEORY 

The theory of loss cone structure was originally de- 
veloped to model the rates of tidal disruptions of stars 
by MBHs ({Frank fc Reed 1197^1 . Results reviewed here 
were derived in the context of g lobular cluster-like stel- 
lar system s containing MBHs ifCohn fc Kulsrudl Il978t 
llpserl jT978t iLightman fc Shapiro! Il977[) where the stars 
are tidally disrupted after coming within the tidal radius 
of the MBH: 



where m* and r* are, respectively, the stellar mass 
and the stellar radius, and r\ is a form factor of order 
unity. Recently, the loss cone paradigm has been uti- 
lized to estimate rates of disruptions of stars by M BHs 
residing in galaxy nucl ei (Magorrian & Trcmainel 119991 
iSigurdsson fc ReesllT997t iSver fc Urmerlll999ft . aswell as 
to calculate the rates of scat tering of stars into the cap- 
ture zone of MBH binaries |Yul l2QQ2t ). At the heart of 
these studies lies the assumption that stellar systems, 
such as globular clusters or small galaxies, are old com- 
pared with local two-body relaxation times. This justifies 
the description of these systems in te rms of the time- 
independent Fokker-Planck equation l ) Cohn fc Kulsrudl 
ll9ralRosenbluth. MacDonald fc Juddlll957|) . 

For simplicity we restrict attention to spherical galax- 
ies; there the phase space distribution is a function of 
two orbital integrals — the energy E and the angular mo- 
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mentum J of a star. The Boltzmann equation reads: 



df 
dr 



v r BR 



((Ai?) 2 ) df 
2R OR 



(2) 



(e.g. iCohn fc Kulsrudl 119781 eq. 30), where 

f(E, R, r)d 3 rd 3 v is the probability of finding a star 
in the phase space volume element d 3 rd 3 v, v r is the 
radial velocity of the star, R = J 2 /J 2 (E) is the square 
angular momentum in units of the angular momentum 
of a circular orbit at energy E, and ((Ai?) 2 ) is the 
diffusion coefficient associated with the orbital integral 
R. At no risk of confusion we will variably refer to R as 
"angular momentum." Note that all terms representing 
diffusion in energy have been discarded; this is justified 
because the gradients in E are minuscule compared 
with the gradients in R in the vicinity of the loss-cone 
boundary. The diffusion coefficient scales inversely with 
the number of stars in the galaxy: 



((Ai?) 2 ) oc N- 



(3) 



1 galaxy 

Conditions for the validity of the time- independent, 
orbit-averaged Fokker-Planck equation can be summa- 
rized as follows: 



P{E) « T rc i ax (£) « T a; 



(4) 



where P(E) is a typical orbital period (or the crossing 
time) of a star at energy E, T^iax is the local two-body 
relaxation time, and T agc is the total lifetime of the sys- 
tem. If the first inequality fails, perturbative expansion 
of the Boltzmann equation is no longer valid, large de- 
flections dominate, and the system is in a gaseous phase. 
If the second inequality fails, the system need not be col- 
lisionally relaxed and can keep evolving; then the time- 
dependent Fokker-Planck equation is necessary. We ar- 
gue in the next section that the second inequality is al- 
most never satisfied in galaxies with MBH binaries. 

Let fbin denote the distance from the binary's center of 
mass within which stellar orbits are strongly perturbed 
by the rotating quadrupole gravitational field of the bi- 
nary. This distance is similar to the semi-major axis 
of the binary, rbi n ~ a- Stars that transgress a region 
around the binary r < r^n are subject to gravitational 
slingshot and are ejected from the system. This imposes 
a boundary condition on equation J2J) of the form: 

f(E,R,r)=0 (r<r bin ). (5) 

Ligh tman fc Shanircl l|1977f) derived an orbit averaged so- 
lution f(E,R) to equation J5J subject to the boundary 
condition expressed by equation (jSJ) to find: 



R 

Rq 



(G) 



f(E, R) ee P-\E) I -f(E, R, r) = Cin 

J V r 

where C is a normalization factor to be determined, and 
i?o (E) is an energy-dependent angular momentum cutoff. 

Equation © assumes that f(E,R,r) vanishes identi- 
cally when R < Rq. In the absence of two body relax- 
ation, i?o = -Ric, where R\ C (E) is the angular momentum 
of a particle with pericenter distance at the very edge of 
the sphere of capture, r_(_E, R\ c ) — r^ ln . In the presence 
of relaxation, some stars will enter and exit the consump- 
tion zone R < R\ c all within one orbital period, thereby 
evading consumption. ICohn fc Kulsrudl l)1978|) carried 



out a boundary-layer analysis to find the following rela- 
tion between Rq and R\ c : 

f>(F\-R /pwl cx P( _ ?)' Q(E)>1 
H °^> - Klc ^ X \exp(-0.1865-0.824^), q(E) < 1 ' 

(7) 

where q(E) is the ratio of the orbital period at energy E 
to the timescale for diffusional refilling of the consump- 
tion zone at this energy: 



Ric{E) 



dr ((Ai?) 2 ) 

v r R^O 2R 



(8) 



When the period is much shorter than the refilling time 
scale (q <C 1), the system is in a "diffusive" regime, the 
loss cone is largely empty and Rq ~ i?i c - When the 
period is much longer than the refilling time scale (q 
1), the system is in a "pinhole" regime, the loss cone is 
largely full and Rq <C i?i c . The energy at which q = 1 is 
called the critical energy. If the loss cone is as narrow as 
in the tidal disruption problem, parts of the phase space 
at the largest energies E > 2a 2 will be in the diffusive 
regime, while the less bound regions will be in the pinhole 
regime. Most of the lost cone flux originates near the 
critical energy. 

The normalization factor C can be expressed in 
terms of the isotropized distribution function f(E) = 

/g 1 f(E, R)dR to obtain: 



C(E) = (R Q -lnR -l) 1 f{E), 



(9) 



implying that the orbit-averaged solution for the time- 
independent Fokker-Planck equation reads (using Rq <SC 
1): 

f^R) = r^^f(E). (10) 
ln(l/i? ) - 1 

In this case, the number flux of stars into the loss cone 
is given by: 



T(E)dE = 4ir 2 J? 



f 



' ln(l/i?o) 



dr ((Ai?) 2 ) 
— hm 

v r R-^o 2R 



-dE, 



1 



(11) 



where all quantities on the right hand side depend on 
energy, and Rq also depends on the geometric size of 
the BH binary, Rp ~ G(Mi + M 2 )a/J 2 . Note that the 
product ((Ai?) 2 )/ is independent of the number of stars 
in the galaxy N ~ Mgaiaxy/m*. The total mass flow 
J T{E)dE, however, scales as A^ 1 . 

3. DOES THE STANDARD THEORY APPLY TO MBH 
BINARIES? 

The orbital separation a of MBHs in a binary is much 
larger than the tidal disruption radius around an isolated 
MBH. Using equation 0J: 



^«10 5 xrr 2 , 



2/3 ( M. 



-1/3 / x 1/3 

Mo 



Re. 



1 pc 



(12) 



Henceforth, M, = Mi + M 2 denotes the total mass in 
MBHs. In the context of a loss cone associated with a 
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binary, most of the galaxy is in the strongly diffusive 
regime and the boundary layer is of negligible thick- 
ness, Ro ~ R\ c . Figure shows the dependence of 
the quantity q on energy for a model described by the 
p ~ r _2 (l + r 2 )~ 2 stellar density profile and scaled to the 
galaxy M32 with 0.1% of the galaxy mass in the central 
MBH. Note that q{E) < 1 when the MBHs first form 
a bound pair and remains <C 1 even when the binary 
separation shrinks to 0.01 pc. Scaling to other galaxies, 
assuming a nuclear density profile of p ~ r~ 2 and a fixed 
ratio of MBH mass to the galaxy mass, we find: 

with a the ID stellar velocity dispersion. Since q(E) ~ 
P(E)/T Ie \ ax (E), the typical deflection of a star due to the 
presence of other stars over one orbital period is small 
compared with the size of the loss cone. Therefore, / 
is nearly uniform around the orbit of the star; the stel- 
lar system satisfies Liouville's the orem; and the orbit- 
avera ged Fokker-Planck equation : l.iuin man fc Sh apiro 
is a valid description of the loss cone dynamics. 
We will rely on this conclusion throughout this work. 

Another aspect of the standard theory fares less well 
in the case of MBH binaries. We have so far discussed 
the time-independent Fokker-Planck equation; strictly 
speaking, however, time-independent solutions do not 
exist as stars that diffuse into the loss cone are moved 
from one orbit to another, or removed from the system. 
The flux of stars into the loss cone given by the stan- 
dard theory is shown in Figure As the stellar orbital 
integrals diffuse, the potential changes accordingly, and 
the isotro pic distribution f(E) adjusts t o the changing 
potential. Mil osavlievic fc MerrittJ l|2001|) demonstrated 
that the change of the overall density profile of the galaxy 
can amount to a destruction of the steep central cusp. 
One may attempt jto account for the changing galaxy pro- 
file by adjusting f(E) as dictated by the isotropic Jeans 
equation, while keeping the formal solution in equation 
HI OB unchanged. Similarly, as the MBH binary hardens, 
the loss cone boundary R\ c decreases. One could try 
scaling the formal solution to accommodate a changing 
Ric- 

This strategy meets a serious obstacle when consider- 
ing the time required for the product of a galaxy merger 
to converge to the quasi-steady state described in the 
previous section. The merger and the subsequent for- 
mation of a hard binary proceed on the local dynamical 
time scale, which can be much shorter than the colli- 
sional relaxation time scale. Therefore the distribution 
function f(E, R) immediately following the formation of 
a hard binary can be far from the solution in equation 
HI 0(1 . The mass in stars near the loss cone in excess of the 
equilibrium solution can be orders of magnitute larger 
than the mass fed into the loss cone over a Hubble time 
assuming the equilibrium solution. Sudden draining of 
the loss cone during formation of the hard binary pro- 
duces steep phase space gradients that are closer to the 
step function: 

f(E,R)^S.f( E ^ . (14) 



Since the collisional transport rate in phase space is pro- 
portional to the gradient of / with respect to R, steep 
gradients imply an enhanced flux into the loss cone and 
an accelerated evolution toward the equilibrium form. 




Fig. 1. — Characteristic quan tities a ssocia ted with loss 
cones in galaxies described by the Uaffd 119831) model, p ~ 
r~ 2 (1 + r 2 ) - 2 . Plots have been scaled to match the galaxy 
M32 with a 3 X 1O 6 M central MBH. The energy is expressed 
in terms of the central ID velocity dispersion a. (a) The func- 
tion q(E). From bottom up, solid curves correspond to loss 
cones with capture radii of 1 pc, 0.1 pc, and 0.01 pc, respec- 
tively. Note that q(E) -C 1 in all cases, indicating that the loss 
cone associated with a binary MBH at the center of a galaxy 
like M32 is always in the diffusive regime, (b) The flux T of 
stars into the loss cone assuming the standard result reviewed 
in §0 the flux is expressed in MBH masses per unit E/2a 2 per 
Myr. (c) The time scale in Myr on which collisional relaxation 
can cause a significant change in the shape of the stellar distri- 
bution function f(E,R), assuming that the all stars have the 
same mass, m» = IMq. Convergence of arbitrary i nitia l con- 
ditions to the standard loss cone solution (equation 1101 takes 
place on this time scale. 
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Figure^ shows that in low mass ellipticals like M32, 
the time needed to achieve equilibrium depends on the 
energy, but still measures well in excess of the Hubble 
time for the relevant range of energies. In intermedate- 
mass and massive galaxies, the time is longer than the 
Hubble time regardless of the orbital energy. In § 0] we 
derive equations describing the evolution of the stellar 
distribution starting from arbitrary post-merger initial 
conditions and quantify the discrepancies with the steady 
state solution. We find that the steady-state solution is 
never a good description of a stellar system containing a 
binary MBH. 

We have so far ignored the fate of the stars that are 
ejected by the binary. Unlike the case of tidal disruption 
of stars by a MBH, a star ejected by a binary MBH via 
the gravitational slingshot remains inside the galaxy and 
its host dark matter halo. Ejected stars may therefore 
return to interact with the binary for a second and subse- 
quent times. Every ejected star can, in principle, interact 
with the binary infinitely many times. Its energy can in- 
crease or decrease at every encounter. This process has 
been ignored in the published Fokker-Planck treatments 
of tidal disruption; nonetheless it warrants consideration 
because of its capacity to shorten the lifetime of a MBH 
binary even after all of the stars inside the loss cone have 
been ejected at least once. We discuss the efficiency of 
reejection in § El 

The standard treatment predicts that the decay time 
of MBH binaries scales with the inverse stellar mass 
a/a ~ m" 1 ~ N, where N is the number of stars in 
the galaxy. Evidence that self-consistent TV-body simu- 
lations reproduce this trend has so far been equivocal. 
Simulations wit h M./m* < 10 3 (N < 10 5 ), such as 
those reported in Milosavlicvic & Mcrritt (2001), appear 
to have taken place in the "pinhole" regime, where the 
loss cone is almost full and the TV-dependence is weak. In 
simulations where the MBH binary was allowed to wan- 
der in space, its Brownian motion was a likely source of 
additional, spurious diffusion of stars in the phase space. 
Simulations with M./to* ~ 10 4 and a binary that does 
not wander are able to recover the trends characteristic 
of the diffusive regime. We discuss iV-body simulations 
of MBH binary dynamics in § and § [7| 

4. TIME-DEPENDENT EVOLUTION OF THE LOSS CONE 

4.1. The Loss Cone as an Initial Value Problem 

In this section we study the diffusion of stars into the 
loss cone of a MBH binary following a merger of two 
galaxies. We assume that stars are removed from the 
final galaxy the first time they are ejected by the MBH 
binay; the continued role of these stars is addressed in 
§ El The loss cones of galaxies harboring MBHs are in 
the diffusive regime. A large fraction of stars scattering 
into the geometric loss cone are kicked out by the binary 
in one orbital period. This suggests that R = R\ C (E) 
can be thought of as an absorbing boundary condition 
where the distribution function vanishes. Since R\ C (E) 
has no explicit dependence on the orbital phase, one can 
seek solutions of the Fokker-Planck equation that have 
no explicit dependence on r. The only remaining de- 
pendence on r enters through the diffusion coefficients 
in the collision term of the Boltzmann equation and can 
be integrated out. As a result, one obtains the "orbit 
averaged" Fokker-Planck equation (Lightman & Shapiro 



1977). 

Instead of the distribution function /, we work in terms 
of the number density of stars in the (E, R) plane, which 
is related to the distribution function via: 

N(E, R, t)dEdR = 4tt 2 P(E, R)J c (E) 2 f(E, R, t)dEdR. 

(15) 

The time-dependent Fokker-Planck equation including 
terms that describe diffusion in the angular momentum 
direction reads: 



9N 8 ,, An , lft 1 d 2 



[((AR) 2 )N] . (16) 



Using a relation ijBinnev &: Lacevlll988j) between the the 
first- and the second-order diffusion coefficients, 



equation (|16fl can be simplified to obtain: 



dN 
~dt 



ld_ 

2dR 



Next we expand the diffusion coefficient ((AR) Q 
limit R 0: 



\((AR) 2 ) = 



lim 

R— »0 



<(Ai?) 2 



2R 



R + 0(R 2 



(17) 



(18) 



in the 



(19) 



where the coefficient in brackets depends on energy 
and radius. We then average over one orbital period 
P^ 1 § dr/v r so that the Fokker-Planck equation (R <C 1) 
becomes: 



dN 
~dt 



d 



R 



dN 



OR V dR 



where 



1 /*,. ((AR 2 

u(E) = — <f> — lim ^ '-!- 

PV ' AP J v r R^o 2R 



(20) 



(21) 



Note that this definition of p ,(E) differs from that in 
Magorrian & Tremaina l)1999j) by a factor of 1/4. We 
finally carry out a change of variables R = j 2 to rewrite 
the equation as: 



dN 



fj, d ( .dN 



j dj V dj 



(22) 



This is the heat equation in cylindrical coordinates with 
radial variable j and diffusivity [i. At every energy E, 
boundary conditions need to be specified at < j\ c < 1 
and j = 1, where ji c = y/Ri c - The boundary condition 
at j = 1 is of the Neumann type: 



dN 



= 0. 



(23) 



The boundary condition at j — j\ c is a perfectly absorb- 
ing, Dirichlet boundary condition: 



N(E,j)=Q (j<j ic ) 



(24) 



The geometrical meaning of the boundary conditions is 
illustrated in Figure [21 
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J = l 




Fig. 2. — Schematic representation of the loss cone. The 
orbital population function when the hard MBH binary forms, 
Af(E,j, t = 0), is specified across the shaded surface (cylindrical 
symmetry is assumed). For t > 0, the evolution of Af is identical 
to the diffusion of heat into the sink region j = j\ c . As the MBH 
shrinks, ji c moves inward. 

The general solution of equation (|22|> subject to bound- 
ary conditions l|23|) and l|24|) can be obtained by means 
of Fourier-Bessel synthesis (e.s. lOzisiki ri98C0 : 

J\f(E j t) = — V [PmMPmjlc)} 2 

xA(/3 m ;j)e-^* 

x / j'A(P m ;j')M(Ej',0)dj', (25) 



where J n and Y n (n — 0, 1) are Bessel functions of the 
first and second kind, A(x; y) is a combination of the 
Bessel functions defined via: 

A(x;y) ee J (xy)Y 1 (x) - J x {x)Y {xy), (26) 

while the /3 m are consecutive solutions of the equation: 

A(p m ;j lc )=0. (27) 

We also evaluate the time-dependent version of the loss 
cone flux: 



T(E,t)dE: 



d 
dt 



2 j Af(E,j,t)jdj 

jlc 



^ 2 E 



=i [MPra]l C )\ 2 ~ [JliPm)] 2 

xB(/3 m ;j lc )e^^" t 

x / j'A(0 m ;j')Af(E,j , ,Q)dj'. 
Jjic 



(28) 



where B(x; y) is another combination of the Bessel func- 
tions 



B{x;y) ee J 1 (xy)Y i (x) - J^Y^xy). 



(29) 



4.2. History of the Loss Cone Profile 

In Figure |2t we show the evolution of Af(E, R, t) at 
one, arbitrarily chosen E, assuming that the loss cone is 



empty within R\ c and that is a constant function of 
R outside of Ri c at the beginning (t = 0). Here we have 
assumed that the loss cone boundary is static, R\ c ~ 0.02. 
For comparison, we also show the equilibrium solution of 
equation (|10|l which can be expressed: 



A/" C qui(S, R, t) 



HVRic) - 1 



Af(E,R,t)dR. (30) 



The phase-space gradients dAT/dR decay rapidly at first 
and then more gradually as they approach the equilib- 
rium solution; the former converges to, but never be- 
comes equal to the latter. It is evident that the total 
population J MdR incurs a decrement of order unity be- 
fore it has had time to reach the state of collisional equi- 
librium. A decrease in Af implies a decrease in the spatial 
density gradients in the galaxy, or "cusp destruction." 
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Fig. 3. — (a) Slices of the density Af(E, Ft, t) at one, arbitrary 
E, recorded, from left to right, at 10°, 10 1 , 10 2 , 10 3 , and 10 4 
Myr (solid curve). Initially, Xf(E,R,t) = for R < R lc and 
Af(E,R,t) = const fo r R > i?i c . We also show the equilibrium 
solution of equation I1UI (dot-dashed curve), (b) The total 
number of stars consumed by the loss cone as a function of time 
(solid curve). The scale has been set to galaxy M32 modeled as 
a Jaffe model Uaffeil983l) with M. = 3 X 10 s Mq and an initial 
separation between the MBHs of 0.1 pc. 

In Figure we present the total mass consumed 
by the binary, which is given by equation (|32f) . The 
time-dependent and the equilibrium masses are shown 
for comparison. The equilibrium solution (dot-dashed 
curve) yields a loss cone flux smaller by at least a fac- 
tor of two within a Hubble time (Fig. |3t>) - The static 
loss cone boundary a ppro ximation employed here is valid 
when Mi os t <C M, f§ !4.3[l . or below the dashed line in the 
figure; we consider the general case in § 14.41 Our model 
here is that of a galaxy with a steep density profile such as 
M32. The enhancement over the equilibrium fluxes is the 
strongest in the galaxies with shallower central profiles 
or "cores;" there, however, Mi os t remains smaller than 
M. during a Hubble time assuming perfect sphericity of 
the galactic potential. 

The following remarks illustrate how the time depen- 
dence of the loss cone profile complicates the inferences 
that can be made about MBH binaries in real galax- 
ies. Given p(r), the Jeans equation can be solved for 
the isotropic distribution function f{E) and the orbital 
population Af(E) = j Af(E,R)dR. However from equa- 
tion (|20() , the flux of stars scattering into the loss cone 
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is proportional to the phase space gradient in R at the 
boundary of the loss region, or: 



T(E, t) = ifiRu 



dN 
~dR 



(31) 



If the nucleus is not collisionally relaxed, the angular 
momentum dependence of N(E, R) is unknown. Since 
the gradient dN/dR evolves in time (e.g. equation l2"%|) . 
the loss cone flux is likely to be a strong function of the 
age of the galaxy since the most recent merger event, as 
well as of other details like the mass ratio of the merger. 

Attempts to estimate the separati ons of MB H binaries 
in presently observed galaxies (e.g. lYul 120021) are there- 
fore problematic. As a first step, we would need to relate 
the total mass in stars ejected by the MBH binary (fol- 
lowing the initial emptying of the loss cone) to the evo- 
lution of its semi-ma jor axis. The most commonly used 
form of this relation l)Quinlanlll996j) . derived via scatter- 
ing experiments and confirmed i n TV-body simulations in 
the c ase of equal- mass binaries ( Mil osavlievic fc Merrittl 
l200l|) . reads: 



M\ ost (t) =171* 



< JM. In 



F{E,t')dEdt' 
a(0) 

a{ty 



(32) 



where, again, M. = M1+M2 is the binary mass, and J is 
a numerical factor approximatel y equal to one when the 
binary is hard and M\ ~ M 2 . iMilosavlievic fc Merritti 
(2001) found that J ss 0.5 correctly described the bi- 
nary decay in their A-body simulations; in a general case 
J is a function both of the hardness of the binary and 
the MBH mass ratio Mi/M 2 (but see §E2J. Therefore 
the number of e-foldings in the binary separation is pro- 
portional to the mass ejected. Since the evolution of the 
binary separation depends on the history of the loss cone 
flux from the binary formation onward, the assumption 
that the equilibrium flux obtained from the presently ob- 
serve d p(r) was the same at all times in the past (jYul 
2002) is too conservative. Not only would the phase- 
space gradients dN/dR have been steeper in the past, 
but the density profile would have been steeper as well, 
providing significantly enhanced fluxes early in the life 
of the binary. Present-day MBH binar y separati ons cal- 
culated assuming steady-state galaxies l)Yull2002ft should 
therefore be interpreted as upper limits. 

4.3. Approximations 

Several approximations have been made in deriving the 
solution in equations 125|) and l|28l) . These approxima- 
tions relate to the assumption that the isotropized dis- 
tribution function f(E) changes very little over the time 
interval t G [0,i max ] when the solution is expected to be 
valid. We summarize the approximations as follows. 

1. The local crossing time is much shorter than the 
local relaxation time, P{E) <C T le \ ax (E). This is gen- 
erally required of the perturbative expansion on which 
the derivation of the Fokker-Planck is based. Since we 
have calculated the time dependence explicitly, it is no 
longer necessary that the system be old compared with 
the relaxation time. It is however imperative that the 
physical system modeled by equation i|22|) have reached 



dynamical equilibrium prior to the time the model can 
be considered applicable. 

2. The shallowing of the overall potential due to the 
ejection of stars from the loss cone is small. This guar- 
antees that the diffusivity n(E) also stays constant. In 
the tidal disruption problem this is always the case, as 
the consumed stellar mass is at most a small fraction 
of the MBH mass. In the massive binary problem, the 
potential is never constant, especially following mergers 
of galaxies with steep cusps when the binaries impart 
significant damage to the merged cusp. Indeed, if the 
ejected mass becomes comparable to the MBH mass — 
which is a likely outcome of mergers of galaxies with steep 
density cusps — the potential near the MBH will incur a 
relative change of order unity. This notwithstanding, if 
the potential changes slowly enough, adiabatic invari- 
ance suggests that the solution of equation could 
be very close to the true solution modulo an appropriate 
reparametrization of the energy dependence. 

3. The geometric loss cone boundary R\ C (E) = j\ c 2 {E) 
changes little in time. Again, this is true for the tidal 
disruption loss cone, but almost never true for the bi- 
nary loss cone. As the binary siphons energy into stars 
that are captured and then ejected, the semi-major axis 
a of the binary decreases. Since R\ c ~ a 2 , j\ c shrinks 
proportionally to the binary separation. The solution 
in equation (|25(l will be valid only so long as Aa <C a. 
The domain of validity of the static loss cone boundary 
approximation can be estimated by considering M\ ost (t). 
Recall that the loss cone boundary scales with the bi- 
nary separation, R\ c ~ a. Therefore the approximation 
in which the geometric loss cone boundary is static is 
valid only if Mi ost (i) < JM,. 

4.4. Evolution of the Semi-Major Axis 

A self-consistent solution accounting for the evolution 
of the loss-cone boundary in response to the binary's or- 
bital decay can be constructed by iterating over the ap- 
propriately chosen time steps. The steps are chosen such 
that the binary separation changes little within each step, 
At < JM,/m* J T[E, t)dE. At the end of each time 
step, one corrects the binary semi-major axis a for the 
energy that the binary has exchanged with the ejected 
stars as dictated by equation Ij32jl . One then recalcu- 
lates the Fourier-Bessel coefficients in equation Ij25[) us- 
ing the new, advanced N(E, j, t + Ati) and the corrected 
loss cone boundary j\ c (E,t + Ati) and thus extends the 
time-dependent solution beyond t + Ati. Iteration of this 
kind overcomes the requirement for a static background 
and presents a computationally feasible numerical proce- 
dure for calculating loss cone dynamics over cosmological 
times. 

We carried out this procedure in a model of the galaxy 
M32 and present the results in Figure 0Ji-c. The char- 
acter of a(t) does not vary with the choice of the ini- 
tial separation et(0) as long as the initial width of loss 
cone is consistent with the separation chosen; here we 
have set the geometric loss cone boundary to equal 
Ric{E) = GM.a/ 'Jc(E). The initial separation used in 
the middle panel of Figure a(0) = 0.1 pc, is the clos- 
est to what we expect to find in the merger of two equal 
galaxies with ~ 10 6 M Q MBHs. Note that while \da/dt\ 
is bounded in the equilibrium solution (dashed line), in 
the time-dependent solution (solid line) \da/dt\ is subtan- 
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tially larger for small t and diverges in the limit t — > 
due to the very large initial gradients dN /dR. 
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Fig. 4. — (a)-(c) Evolution of the semi-major axis in a singu- 
lar isothermal sphere scaled to galaxy M32 with initial values 
of, from top to bottom, 1 pc, 0.1 pc, and 0.01 pc. Solid line 
is the evolution predicted by equation 1321 ; dashed line is the 

:§E». 



prediction of the equilibrium theory (§ 



Scaling to other 



galaxies with p ~ r stellar density cusps can be achieved by 
linear re-parametrization of the time, (d)-(f) Enhancement of 
the exact solution over the equilibrium solution (triangles). The 
quantity © is defined in the text. Power-law fits to the ©(i) in 
the first 1 Gyr are also plotted (solid lines). 

We can express the difference using the quantity (de- 
noting the equilibrium solution with a oqu i): 

e(t) = °f~ 0( *L (33) 

a[0) - a cqui (t) 

In Figure 0Ji-f we show the evolution of O(i); the en- 
hancement is well described by: 

3xlO 6 M m* t 



m~ao-is)x y M% M0lOMyr 

This relation is a good fit for t < idostr where t<jestr is 
the time scale on which the removal of stars from the 
stellar cusp (the cusp destruction) starts affecting the 
pool of stars available for diffusion into the loss cone. 
The time-scale can be identified with the diffusion of the 
stellar mass of about JM, ~ Q.5M. into the loss cone, 
which in turn corresponds to the change of the binary 
separation by the factor e -1 / 2 ~ 0.6. From Figure0^,-b, 
a/a(Q) — 0.6 at approximately at 1 Gyr for a 3 x 10 6 M Q 
MBH. Scaling to an arbitrary binary mass gives: 



— x 1 Gyr. 



(35) 



3 x 10 6 M Q 

The difference between the two solutions decreases with 
time, to roughly a factor of two after a Hubble time. 



There is, however, a plausible scenario in which the 
mechanism discussed in this section could lead to a much 
greater enhancement in the binary's mean rate of decay. 
Galactic nuclei may episodically return to a strongly non- 
equilibrium state due to events like the capture of a dwarf 
galaxy, infall of a giant molecular cloud, or star forma- 
tion from ambient gas. A rejuvenated loss cone will once 
again result in strong phase-space gradients as the bi- 
nary rapidly ejects stars with R < Ri c . Then the the 
non-equilibrium enhancement averaged over the entire 
lifetime of the binary would be higher than suggested by 
Figure^] For example, if the episodic replenishment of a 
MBH binary loss cone in a galaxy like M32 occurs every 
10, 100, or 1,000 Myrs, equation (J31J implies that the de- 
cay of the binary in that galaxy will occur, respectively 
O w 10, 5, or 3 times faster than what the equilibrium 
theory would have predicted. 

In conclusion: the equilibrium solution in equation 
111 Oil — originally developed to describe the much smaller 
loss cones representing tidal disruption of stellar mass 
objects by massive ones — is never a good approximation 
to the structure of binary MBH nuclei. Instead, an epoch 
in which the flux of stars into the loss cone is enhanced 
compared to the equilibrium solution, is succeeded by 
one in which the cusp is significantly modified due to 
ejection of stars. 

5. SECONDARY SLINGSHOT 

5.1. General Considerations 

Until now we have ignored the possibility that the stars 
ejected by the binary once are ejected again as they re- 
turn to the nucleus on radial orbits. The stationary so- 
lution for the loss cone dynamics (§ 0) and the time- 
dependent solution (§ 0J are based on the "sinkhole" 
paradigm modelled after the tidal disruption of stars by 
a MBH. In the tidal disruption picture stars are destroyed 
as soon as they transgress the tidal disruption zone of the 
MBH. This is not true for MBH binaries. The capture 
cross section of binaries is larger than that that of sin- 
gle MBHs by factors of 10 5 — 10 11 and only a negligible 
number of stars are disrupted; the vast majority simply 
receive a kick (AE, AR) that transports them to another 
orbit. Following an event of slingshot, the final angular 
momentum is comparable to the binary's orbital angu- 
lar momentum, J ~ \ZGAI t a. If the binary orbit decays 
slower than the orbital period, (dlna/dt)^ 1 P, most 
stars inside the loss cone remain inside the loss cone, en- 
counter the binary at their next pericenter passage, and 
are "reejected." If the decrement in the binary separa- 
tion is substantial, some stars that are ejected very near 
the loss cone boundary finish just outside boundary as 
the boundary shifts inward. These stars are lost at a 
rate proportional to the rate of decrease of the loss cone 
angular momentum R\ c cx a. 

In this section we explore the possibility that the 
reejection can prolong the decay of a MBH binary beyond 
the stalling point when the loss cone would have been de- 
pleted in the sinkhole paradigm. To simplify the forth- 
coming derivation, we consider an orbit-averaged system 
in which probability for slingshot of star of energy E lo- 
cated inside the loss cone is uniform in time and totals to 
unity within one orbital period. Some stars can become 
quasi-permanently captured by the MBH binary; in re- 
ality such captures are possible but statistically negiligi- 
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ble. Furthermore, we assume that the orbital population 
function Af(E, R 7 t) is constant in R inside the loss cone. 
We define the loss-cone orbital population function: 

/■flic 

Af(E,t) = Af(E,R,t)dR. (36) 



To derive the time-dependence of J\f(E,t), we consider 
the evolution over a small interval At. The evolution 
has to account for the redistribution of stars due to the 
kicks imparted by the MBH binary, as well as for the 
gradual influx of stars due to the two-body relaxation. 
Therefore we have: 

Af(E, t + At) = J Af(E', t)C(E', E, At)dE' + F{E, t)At. 

(37) 

Here, Q{E' , E, At) is the transition probability between 
the energy E' at time t and the energy E at time t + At, 
and T[E, t) is the relaxation flux of stars as derived in 
§ |2 and § 0] Note that both depend implicitly on the 
binary's semi-major axis. The transition probability can 
be understood as a combination of the kicks from E into 
some other energy, the kicks from another energy E' into 
E, and the attrition related to the decrease in the loss 
cone size R\ c : 



C(E',E ) At) = 



1 - 



At 



P(E') 
At 



S(E-E') 



P(E') 

R\ c {t 



(i(E',E) 
At) 



«.w "J ' (38) 

Here, (i(E', E) is the probability that a star ejected from 
energy E' will end up on an orbit with energy E. The 
probability is normalized to unity J (i(E' , E)dE' = 1. 

Equation l|38|) is written under the assumption that the 
orbital phases of stars at a given energy are randomly 
distributed. In reality, this may not be the case as stars 
retain the memory of their initial orbital phases. Vari- 
ance in the ejection velocities will generate dispersion in 
the distribution of stellar phases, which will destroy the 
phase coherence of the stars after they are ejected a few 
times. We do not believe that an initial non-uniform dis- 
tribution of the phases will have a significant effect on 
the decay of black hole binaries. 

Passing to the infinitesimal limit in equation (|37|l and 
using d/dt\\\R\ c (t) — d/dt\na(t) we obtain an equation 
of motion for the orbital population of the loss cone: 

d M(E,t)_ JV(E,t) , [«{E\t) <i[E ^ E)dE , 



dt 



P(E) 
-Af(E, t) 



P(E') 



dlna(i) 
dt 



T{E,t). 



(39) 



We proceed to derive an equation describing the evolu- 
tion of the total number of stars inside the loss cone: 
dM(t) _ 



dt 



j t J Af(E,t)dE 



=Af(t) 



dlna(t) 
Jt 



J F{E,t)dE, 



(40) 



and of the total energy budget of the loss cone: 
d£(t) _ d 
dt = ~dt 



Af(E, t)EdE 



j^0\-E + j^E,Ei)E>dE> 



P(E) 

, d In ait) 
Af(E,t) — j^EdE 



dE 



J T{E 1 t)EdE. 



(41) 



The physical meaning of the factor in brackets in the last 
row of equation 141fl is immediately apparent — it is the 
average energy change that a particle originally at en- 
ergy E experiences as a consequence of the gravitational 
slingshot: 



(AE) = -E - 



J Ci(E,E')E'dE'. 



(42) 



The energy gained in the first term of equation 1411) has 
to be compensated by the change in the binary's binding 
energy 



d (GM X M 2 
~dt V 2a 



P(E) 



(AE)dE. (43) 



The hardening of a MBH binary coupled to an evolving 
stellar population inside the loss cone is described by 
equations Ip^l and l(i3^l. These equations are rendered 
in a form conducive to numerical integration. 

5.2. Example: The Singular Isothermal Sphere 

To obtain a sense of how the secondary slingshot con- 
tributes to the decay of the binary, we consider the 
following example: the galaxy is a singular isothermal 
sphere (SIS) with the density law: 

M 

PM = 7ZT^> (44) 



47T7'o r 2 



and the potential: 



$(r) = In 



(45) 



where M and ro are constants. Note that the potential 
due to this configuration diverges at r — * oo; we have 
been free to fix $(ro) = 0. 

The SIS is probably a better approximation to the po- 
tential of a generic early-type galaxy than it may seem. 
Indeed, while the luminous matter (the stars) in a galaxy 
follows the SIS profile only within a fraction of the effec- 
tive radius, the combined stellar and CDM contributions 
to the potential appear to approximate the SIS over and 
beyond the extent of the luminous distribution. For sim- 
plicity we also assume that relaxation is absent in the 
model: T(E, t) = 0. The radial period inside the SIS 
potential is: 



-E/2a 2 



(46) 



where a = y / GM/2r is the velocity dispersion outside 
the radius of influence of the MBH. Substituting this into 
equation (|4^|l gives: 

(47) 
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Naively it may appear that the integral in equation l|47|l 
diverges. This not the case because large (very bound) 
energies are forbidden inside the loss cone, Af(E, t) = 
for E > E max , since the coexistence of a star with large 
binding energy and a binary is very improbable. The 
distribution of stars inside the loss cone peaks at the 
energy corresponding to the slingshot ejection velocity 
(see Figure [SJ: 



E, 



eject 



eject 



1 G(Mi 

2 c 



Af 2 ) 



(48) 



where $ eject is the fiducial potential energy of a star at 
the edge of the MBH binary's sphere of influence and 
v /G(Mi + M 2 )/a is the orbital velocity the MBHs. To 
be able to solve equation l|47|l analytically, we pretend 
that all stars inside the loss cone have the same energy 
equal to -E G ject at time t: 



Af(E,t) -f Af(t)5(E - ^ject). 



(49) 



Therefore: 



d / GM\M 2 
~dt V 2a 



P(0) 



Af(t)(AE) } 



t/2* 2 



(50) 

Equation (|50|l can be simplified further after analyzing 
the evolution of the product Af(t) (AE) Bcjoct . Equation 
(|40|) implies that in the absence of relaxation, the total 
number of stars inside the loss cone decays proportion- 
ally to the binary separation, Af(t) cx a(t). We further 
argue that (AE)E ej<:at exhibits scaling with a^ 1 (t), by 
considering two extreme cases: 




E/(2a 2 ) 



E/(2a 2 ) 



Fig. 5. — Left panel, evolution of the distribution of stars 
initally inside the loss cone M(E,t) in a simulation with V ~ 
18, 000 stars inside the loss cone at t = 0. The distributions, 
from right to left, were recorded at logarithmically progress- 
ing times t = (0,1,2, 4, 8, 16, 32, 64, 128, 256) x 4 x 10" 4 P(0), 



where P(0) 



rV2 r 



is the radial crossing time at ro. 



The binary separation decayed by the factor of ~ 4.0 during 
the same interval. The MBH binary mass M, is 10 — 4 of the 
SIS mass M, and the initial binary separation is O.lGAf.cr -2 , 
consistent with the separatio ns found in merger simulations 
( Milosavlicvic & Merritt 2001). The smooth curves were gen- 
erated from iV-body data via a maximum penalized likelihood 
technique. Peak of the distribution travels from larger to 
smaller binding ener gies following the increase in the energy 
of ejection (equation 1481 . Note the late formation of the sec- 
ondary hump at E/(2<j 2 ) ~ 3.0 - 8.0, levelling at V ~ 1,400; 
this hump consists of stars that had originally been inside the 
loss cone, but then finished just outside, out of reach to the 
MBH binary. Right panel shows V = 10 3 cW(-E, t)/d hit at 
t = (1, 2, 4, 8, 16, 32, 64, 128) X 4 X 10~ 4 P(0). The zero-intercept 
indicates the energy at which the average exchange of energy 
with the MBH binary vanishes, (AE) = 0. 

1. Just after two MBH form a hard binary, the stars 
in the inner stellar cusp inside the loss cone interact with 



the MBH binary for the first time and receive a kick com- 
parable to, if not larger than, than their kinetic energy: 



AE > 



AE | 



1G(M 1 +M 2 ) 
2 a 



(51) 



Thus (AE) oc a _1 (t) is indeed a good approximation. 

2. When the binary is old and evolving slowly, one 
may wonder whether a steady state ((AE) — 0) can be 
achieved, in which the stars in the loss cone exchange 
no net energy with the binary. 2 A steady state might 
be expected when the majority of stars inside the loss 
cone have been ejected to large energies at which the 
energy exchange with the binary is very inefficient. As- 
suming that we have indeed succeeded preparing such 
a steady state, we recognize that the second moment 
((AE) 2 ) must be finite. Defining SE = ({AE) 2 ) 1 / 2 > 0, 
we consider the radial periods of stars ejected slower than 
the average (E\ = E + SE) and faster than the average 
(E2 = E — 5E). Since the former return to re-encounter 
the binary sooner than the latter, P(E±) < P(Ei), at 
time t + P(E) the distribution of stars impinging on the 
binary will have acquired a finite first moment: 



E -> E+ (SE)' 



dlnP- 



dE 



E 



(set 

2a 2 



(52) 



Thus the average effective energy of stars as experienced 
by the binary is shifted by (SE) 2 /2a 2 , and the binary will 
respond with a decrease of the semi-major axis by Aa = 
—a 2 (SE) 2 Afm*/GMiM 2 cr 2 . This in turn will generate a 
shift in -E^ect by the amount: 



E, 



eject " 



■E, 



eject 



E, 



G(M 1 + M 2 
2a~ 2 
(SE) 2 



Mr 



eject 



2a 2 



(53) 



where /i = M\Mil(M\ + M 2 ) is the reduced mass. Fi- 
nally, the true distribution of stellar energies shifts in 
response to the incremental binary hardening by the 
amount (AE) = AE C j CCt . Recall that Af oc a. From 
the scaling of the restricted three-body problem we know 
that SE ~ G(Mi + M 2 )/2a. Therefore, again: 



(AE) 



JV(Q)m, a(t) G 2 (M 1 +M 2 ) 2 
fi a(0) 8a 2 a 2 (t) 



(xa^fi). (54) 



We have shown that when the binary just becomes 
hard and evolves rapidly, as well as when it is very hard 
and evolves very gradually, the product Af(t)(AE) cx 
a l a~ l cx a is approximately constant (although the con- 
stant of proportionality could be different in the above 
two regimes). We proceed to integrate equation l|50|l 
where we substitute -E C ject from equation l|48|) : 



1 



1 



Aa 2 



a(0) G(Mi + M 2 ) 



In 



1 



iM(AE) t 
2fia 2 



P(E 

for t > 0, where a(0) is the initial separation and Eq is 
the energy of ejected stars at the outset: 



Eq — E, 



eject \t—o 



eject 



G(Ah + M 2 ) 
2a(0) 



(56) 



We thank C. Pryor for raising this question. 
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Equation l|55[> describes the evolution of a MBH bi- 
nary in the SIS potential ignoring relaxation. In the 
secondary slingshot mechanism, the inverse semi-major 
axis and, hence, the binding energy, are logarithmic func- 
tions time. This result has been derived assuming that 
relaxation-driven loss-cone diffusion was negligible. The 
logarithmic behaviour can be contrasted with the effect 
of diffusion in which the binding energy grows linearly. 
In a realistic galaxy both mechanisms are at work; the 
reejection amplifies the effect of diffusion on the harden- 
ing rate of the MBH binary. 

It is encouraging that the logarithmic behavior seems 
to hold in the simulations where the interstellar interac- 
tions have been replaced by a smooth potential to prevent 
relaxation. Figure shows the evolution of the inverse 
semi-major axis in a simulation with J\f ~ 18, 000 stars 
in the loss cone at t = 0. The observed rate of decay 
da^ 1 /dint « 1.7 x 10 4 comes close to the prediction of 
equation (ESJ: 4cr 2 /G(Mi +M 2 ) = 2x 10 4 . 
4 




t/P(0) 

Fig. 6. — Decay of the inverse semi-major axis l/a(t) in an TV- 
body simulation where the stellar potential has been replaced 
by a smooth component, thereby preventing relaxation (solid 
line). The binary separation does not stall, but continues to 
decay because of the seconary slingsho t (§ 151 . The galaxy is 
a singular isothermal sphere (equation 1441 *1 with M = ro = 1. 



The time is expressed in the units of P(0) 



1-1/2, 



which 



is the radial crossing time at r p. T o compare this curve with the 
model presented in equation 1551 . the hard- binary separation 
a(0) has to be chosen; it will generally be different from a at 
the beginning of the simulation. We show prediction of the 
model with a(0) = 5 X 10" 6 and m*N '(AE)/2fj,a 2 = 2 X 10 5 
(dashed line). 



5.3. Reejection in Galaxies 

An upper limit on the effectiveness of the mechanism 
presented here can be estimated as follows. In a merger of 
two galaxies, the MBHs form a hard binary when a(0) ~ 
Qhard = G^/4er 2 . The total mass of stars inside the loss 
cone just after the hard binary forms is m»A/"(0) ~ 10.0/i. 
We can also assume that Mi = M 2 and (AE) ~ few x 
2a 2 . In a Hubble time, a -1 will have increased by the 
factor: 



a(0) 



a(^Hubble) 



1 + 0.25 In 



few x 10 x 



^Hubble 

P(Eo) 



(57) 



For P(Eq) = 10 3 ' 5 ' 7 years, we obtain a(0)/a(tHubbic) ~ 5, 
4, and 3, respectively. 



This conclusion may be optimistic. In a circular binary, 
orbital velocities of the MBHs can be written as: 



Vi = 2a 



Ohard M 2 



Vo 



a Mi ' " V a M 2 

In the isothermal sphere potential, the stars ejected with 
velocities Vt < v < V2 will leave the binary's sphere of 
influence = GM m /2a 2 and travel to a radius r max : 



exp 



Mg «ha 
Mi a 



< 



< 



exp 



Mi dhaxd 
M2 a 



(59) 



Realistically, reejection becomes ineffective if r max is too 
large: first because the isothermal sphere potential may 
not extend indefinitely; and second because a star with 
large r max is easily perturbed from its nearly-radial orbit 
on the way in or out. If we suppose that r max ~ 100rh 
is an effective upper limit, equation l|5*9l suggests that 
flmin ~ 0.2ohard is an effective lower limit for reejection 
when Mi w M 2 . 

However, we note that many stars are ejected with 
velocities much less than the binary's orbital velocity 
(the velocity change varies between and few x V 2 , and 
can also be negative, depending on the star's orbital pa- 
rameters and the binary's phase). This is particularly 
true when M 2 -C Mi; the average velocity change of a 
star encountering the binary is Av <C V 2 , although every 
star inside the loss cone can be pumped to an ejection 
at ~ V2 after some number of encounters with the bi- 
nary. If Av ~ Vi, which is statistically favored when 
M 2 <C Mi, the lower limit on the semi-major axis is 
weakened, a min <C 0.2ai iar d- We observe that the transfer 
of energy from the binary to the stars is an incremental, 
stochastic process involving multiple encounters; the for- 
malism developed in this section is applicable whenever 
the form of the galaxy's gravitational potential permits 
repeated stellar encounters with the MBH binary. 

5.4. A New Mass-Ejection Formula 

The reejection paradigm suggests a new and more gen- 
eral expression for the relation between the ejected mass 
and the shrinkage factor of a binary. The standard ex- 
pression, equation (|32f> . was motivated by scattering ex- 
periments in which stars are assumed lost if they exit 
the binary's sphere of influence with enough velocity to 
escape the binary. In the case of a binary embedded in 
a galactic potential, the critical quantity is the energy 
gained by a star between the time it enters and exits the 
loss cone. Equating this with the change in the binary's 
absolute binding energy, we find: 

G(M 1+ M 2 ) M ( 1 1 s _ 

— JWiost Renter - -^cxitj- 



itial 



(60) 

In an SIS potential, the greatest contribution to the loss 
cone comes from stars near r^. Defining {E cntcr — E cx i t ) = 
A$ and passing to the infinitesimal limit ainitiai - * 
ag na i = a in equation l|60l) . we obtain: 

1 dMiost 1 , . 

(Mi + M 2 ) dlno-i ~ (A$/2 ( 7 2 )( a /a hard ) 1 ' 

Comparing this relation to equation l|32|l we identify the 
effective value of the mass-ejection parameter J: 

J ~ (A$/2a 2 )(a/ ahard )- (62) 
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To shrink the binary by one e-folding, a stellar mass of 
JM, must be transported from an energy marginally 
bound to the black hole, to the galactic escape velocity. 

5.5. Reejection in Other Geometries 

Reejection would occur differently in nonspherical (ax- 
isymmetric or triaxial) potentials, since angular momen- 
tum would not be conserved and ejected stars could miss 
the binary on their second passage. However there will 
generally exist a subset of orbits defined by a maxi- 
mum pcriccntcr distance that is < a and stars on such 
orbits will interact with the binary once per crossing 
time, just as in the spherical case. As an example, 
consid er centrophilic (chaoti c) orbits in a triaxial nu- 
cleus l|Poon fc Merrittl l2002|L There is one such or- 
bit per energy (in a time-averaged sense), and the cu- 
mulative distribution of pericenter distances for a sin- 
gle orbit is found to be linear in r n up to a maximum 
value, r Piinax (£) l|Merritt fc Poodl2003lK Thus the star 
passes withi n a distance r n max of the center each cross- 
ing time. ijMerritt fc Poonll2003[) evaluated r p<max (E) 
fo r chaotic orbits in the self-consistent triaxial models 
of iPoon fc Merrittl i|2002fl , which have density profiles 
P r* 



and found: 



r P , ma *(E) « 0.3r h e^- E ^ 2 



(63) 



with $ h = $(r h ) and r h = GM,/2a 2 . The probabil- 
ity of a star passing within a during a single traversal 
of the galaxy is ~ a/r p . max for a < r Pimax and one for 
a > ^p,max- It follows that all stars on chaotic orbits 
with E > $h pass inside of ahard on each orbit, and a 
substantial fraction of stars would continue to engage in 
reejection even after the binary had shrunk by a factor 
of several below ahard- Hence reejection in the triaxial 
geometry might occur in roughly the same way we have 
described for spherical systems. The situation would be 
more complex in axisymmetric (nonspherical) galaxies 
but there would always exist some subset of orbits de- 



fined by 



< a particularly if the deviations from 



spherical symmetry were modest. 

6. Af-BODY SIMULATIONS AND REAL GALAXIES 

In this section, we ask whether the long-term evolution 
of MBH binaries seen in TV-body simulations is consistent 
with the theory derived above, and whether the simula- 
tions are capable of reproducing the evolution of binaries 
in real galaxies. 

Does the standard theory (§ |2J reproduce the correct 
time-dependence for the evolution of the binary semi- 
major axis? To solve for a(t), we start from the differen- 
tial form of equation (|32|) : 



dt \aj JM, 



1 dM lost 



[ F(E,a,t)dE. 



dt JM, 

(64) 

Note that both J and T are functions of a. In the 
equilibrium theory feauation lll|) . T oc In -1 (ai/a) where 
a i = J%(E)/GM,e is independent of a. Near the radius 
of influence of the MBH where the pea k loss cone flux 
originates, a\ ~ GM,/2a 2 . Figure 5a of lQuinlanl ljl996f) 
implies that for hard binaries, 



J fa 0.1 x In 



/ 4a h 



arc! 



(65) 



where as before ahard = G/i/Aa 2 . Therefore the evolution 
of the semi-major axis in the standard model has the 
approximate form: 



M 1 

dt \ a 



t/n. 



ln(ai/a) ln(a2/a) 



(66) 



Here, T/N, is a constant proportional to the angular 
momentum diffusion coefficient. We have also explic- 
itly factored out the iV-dependence of this coefficient via 
N, = M,/m*, the number of stars whose mass equals 
that of the MBH. This equation can be easily solved as- 
suming ai ~ a 2 ~ 4ahard, valid for approximately equal- 
mass MBHs (M, = 4/x), yielding: 



°hard _ Ig^ri/AT.) 1 / 3 

a(t) 4 



(67) 



For 1 < Tt/N, < 10, the solution (|67|l is accurately ap- 
proximated by the linear function: 



ahard 



Ft 
N~. 



(1 < Tt/N. < 10). (68) 



For larger values of Tt/N,, we are in the regime where 
the density profile of the galaxy incurs significant damage 
and many assumptions of the standard model do not 
apply. 

In iV-body experiments with V ~ 10 5 and N, < 
10 3 l|Milosavlievic fc Merrit tll2001j) we observed a linear 
growth of l/a(t). However, the binary decay in our sim- 
ulations went beyond the range a(0)/a(t) ~ 1 where the 
assumptions inherent in the standard model are satisfied. 
Additional (unpublished) simulations with N ~ 10 6 also 
produced a linear growth of l/a(t). 

Strikingly, however, these simulations did not yield the 
linear de pendence of l/a(t) on N, predi cted by equation 
P3fl fsee lMilosavlievic fc Merrittl <|200lH . Fig. 8b). Fear- 
ing that the Brownian motion of the MBH binary might 
have contributed to the loss-cone diffusion, we repeated 
the simulations in a symmetric mode that eliminated the 
binary's Brownian motion. 3 The repeated experiments 
did yield a decrease in the hardening rate with V., but 
it was still much weak er than l inear. Similar behavior 
had been reported by iMakinol 119971) whose simulated 
MBHs were more massive relative to the galaxy than 
ours. Makino observed that the decay timescale was sub- 
linear in the number of bodies, da -1 /dt oc V. -1 / 3 , which 
is still strongly at odds with the prediction of equation 

(SSI- 

The explanation for this discrepancy can be found in 
Figure [7| showing the quantity q(E) in these simula- 
tions. Recall that q(E) <C 1 means that the loss cone 
is almost empty ("diffusion regime") while q(E) ^> 1 im- 
plies a nearly- full loss cone ("pi nhole regime ). The fig- 
ure sh ows that the simulations of Mil osavlievic fc Merrittl 
(2001) took place almost entirely in the pinhole regime 
- the loss cone was nearly full at all times. In § 01 we 
argued that the MBH binary loss cones in real galaxies 
are in the diffusive regime throughout the relevant range 
of energies. In the full loss cone regime, the flux of stars 
into the loss cone is given by the number of stars inside 

3 Brownian motion will be identically nil in mergers of equal 
galaxies that are symmetric with respect to reflection (x.v) — » 
(-x,-v). 
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the loss cone boundary assuming no depletion, divided 
by the radial period on which these stars arrive to within 
the binary's sphere of influence, or: 



^(pinhole) ( E } dE 



M(E)R lc {E) 
P{E) 



dE 



(69) 



which is valid for E < E cl - it w 9er 2 (see Fig. 01, 
q(E CT n) = 1. The mass flow into the maw of the binary 
m Jf(^*)(£)d£ is independent of N. and always 
much larger than in the diffusive case. Ignoring reejec- 
tion and noting that R\ c oc a, we find that the semi-major 
axis is predicted to be an exactly linear function of the 
time: 



1 



1 



+ t 



GmJf(E)dE 



(70) 



o(0) ' "Jo JP(E)J 2 (E) ' 

as observed in the AT-body simulations ijMakinol 119971 
Milo savlievic fc Merrittll200H) . 




E/(2a 2 ) 

Fig. 7. — The quantity q(E) in the si mulated merger rem- 
nant A2 of IMilosavlievic L MerritH <200ll) with 131,072 bodies 
before truncation (solid line). Also shown is q{E) in an equiv- 
alent run with 10 6 and no Brownian motion (dot-dashed line). 
The energy is expressed in terms of the central linear velocity 
dispersion a. Note that in the merger remnant q(E) 2> 1 even 
at energies as large as 8<r 2 , implying that the loss cone is in the 
pinhole regime. This is the opposite of what is found in real 
galaxies like M32 that are entirely in the diffusive regime (see 
Fig. and explains why A^-dependence in the simulations was 
weak. 

What value of N is required in order that a nu- 
merical simulation be in the appropriate, diffusive loss- 
cone regime? The scaling of E CI1 t with N can be de- 
duced from equation (|13fl . Assuming a density profile 
p ~ r~ 2 and a potential of the form 2a 2 ln(r/r ) such 
that r = 10 3 GM./2ct 2 we find: 



E, 



crit 



2a 2 In 



7.5 x 10 4 GM./8a 2 
W. a 



(71) 



The transition from a pinhole-like loss cone to a diffu- 
sive loss cone occurs when N m ~ 10 4 — 10 5 . Since a 
typical MBH contains 0.1 % of its host galaxy's mass 
l|Merritt fc Ferraresel l200l . and thus N ~ 10 3 7V., an 
A^-body simulation would have to contain 10 4 ~ 5 x 10 3 = 
10 7 ~ 8 bodies to reproduce the correct, diffusive behav- 
ior of a real galaxy. This requirement is a severe one 
for direct-summation A^-body codes, which rarely ex- 
ceed particle numbers of ~ 10 6 even on parallel hardware 



(e.g. iDorband. Hemsendorf fe MerrittJl2003j) . One route 
might be to couple the special purpose GRAPE hard- 
ware 4 , which is limited to A~ < 10 6 , to algorithms that 
can handle large particle numbers by swapping with a 
fast front end. 

Finally, we can ask how reejection would affect the 
predicted time-dependence of the semi-major axis in N- 
body simulations, assuming a loss cone always full at 
energies E < E cr i t - For simplicity we assume that the 
Brownian motion has been suppressed; the role of Brow- 
nian motion at enhancing binary decay is discussed in 
detail in the next section. Return and reejection of stars 
is a function of the depth of the potential well in which 
the model MBH binary has been placed. The appropri- 
ate mass-ejection coefficient J is the one derived above, 
equation 1(6 2|1 . which takes into account the total energy 
exchange with a star before it finally escapes from the 
galaxy. Note that A<i> — the potential barrier between 
the energy at which stars are fed into the loss cone, and 
the energy at which they escape the galaxy — is a function 
of loss cone entrance energy E. Then: 



at \ a 



f EaHt JT(pinhole)(£) 

J(E) 



M. 

V 

2 \ O-hard 



dE 



(72) 



where, ignoring the implicit dependence of -E C rit on a 
f equation 171(1. I? is a constant: 



V 



G 2 m^ r E ™<- Af(E)A$(E) 
4a 4 J P(E)J 2 (E) 



dE. 



(73) 



The integrand of T> is independent of A. , however -E cr it 
depends weakly on A".. Integration gives: 



flhard 

o(t) 



= Vi + vt. 



(74) 



This result and expression l|55|) were both derived taking 
reejection into account; their difference stems from the 
assumptions made about the rate at which the loss cone 
is being refilled. In equation (|55|l . the refilling is absent 
and the binary decays by reejection only. In equation 
l|74(l . the refilling is at its maximum (as expected in N- 
body simulations) and the reejection serves to amplify 
the effect of refilling on the binary's decay, converting a 
\ogt dependence to a t 1 / 2 dependence. 

7. BROWNIAN MOTION 

Another potentially importa nt source of flux into 
the lo ss cones of simulations in IMilosavlievic fc MerritH 
(2001) was the Brownian motion of the MBH binary. 
Single and binary MBHs alike acquire and maintain a 
random, quasi-periodic drift coming from the discrete 
encounters with stars at the center of the galaxy. The 
drift can be understood in terms of the equipartition 
of kinetic energy between the MBH(s) and the ambi- 
ent stars. The Brownian motion was recently studied 
hv lMerritd ll2001f). iChatteriee. Hernauist fc Loehl (2002) 
and IDorband. Hemsendorf fc MerrittJ l)2003j) . These 
studies concluded that the rms excursions of the MBH's 
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possition and velocity are given by (ignoring factors of 
order unity): 

"' " (75) 



/ 2\ th* 2 



M. 



where r corc is the "core radius" of the stellar distribu- 
tion within which the stellar density approaches a finite, 
central value. In many real galaxies, however, the core 
radius is not well-defined as the density increases near 
the MBH as an approximate power-law, p ~ r~ 7 with 
0<7<2. 

V-body simulations ( Milo savlievic &: Merrittl2001() re- 
veal that the MBH binary remained centered on the stel- 
lar density cusp as it executes the Brownian motion; the 
binary "carries the cusp with it." The cusp endows the 
MBHs with a larger effective dynamical mass. We distin- 
guish the stars closely following the MBH ("satellites"), 
from the stars that orbit within the stellar cusp and come 
close to the MBH during a fraction of their orbital pe- 
riod ("visitors"). The satellites respond to the binary's 
movement and follow the binary as it wanders; their an- 
gular momentum with respect to the binary remains ap- 
proximately constant. The visitors, however, stay em- 
bedded in the static galactic potential and experience 
the binary as moving target. As the binary moves, the 
visitors drift along nearly constant-energy trajectories in 
the E - R plane that can penetrate the loss cone bound- 
ary. The drift can contribute significantly to the loss 
cone refilling flux at energy E if Ai?brown ^ R\c, where 
Ai?brown ~ GM.r\> rown / J 2 (E). 

The meaning of the "core radius" in galaxies with- 
out a core in the density profile can be understood in 
the light of the satellite- visitor dichotomy. The satellites 
drift along with the MBH binary in the background of the 
visitors and other stars not interacting with the binary. 
The satellites, with average velocities larger than the am- 
bient stellar velocity dispersion, form a bound subsystem 
in combination with the MBHs; the subsystem executes 
the Brownian motion in the field of other stars. We can 
therefore think of the background potential in which the 
binary-satellite subsystem moves as the stellar potential 
with the contribution of the satellites subtracted. The 
characteristic radial scale separating the satellites from 
other stars is the the dynamical radius of influence of the 
black holes; therefore the effective value for the "core ra- 



dius" in coreless density profiles such as p < 

GM. 



reads: 



(76) 



Assume the binary executes a small acceleration A, 
which could be the associated with the Brownian mo- 
tion or some other astrophysical process. We think of A 
as a small perturbation to an otherwise orbital integral- 
preserving motion of a test star. In the presence of Brow- 
nian motion, the acceleration can be expressed using the 
quantities presented in equation (|75J) : 



{A 2 



M. r CO r, 



G 2 M? ' 



(77) 



Now let r and J denote, respectively, the position and 
the angular momentum of the star relative to the binary. 
Then: 



dt 



2J 



dt 



(r x v) = -2 J • (r x A). (78) 



The total change in J 2 over one orbital period equals: 

(79) 



AJ 2 = -2 f J ■ (r x A)dt. 
Jo 



Note that AJ 2 (E,3,A) = AJ 2 (E, J,\A\,Cl) is a func- 
tion of the angles Q defining the inclination of the un- 
perturbed orbit relative to the acceleration vector. To 
calculate the second-order Fokker-Planck diffusion coef- 
ficient for the diffusion in the R = J 2 /J 2 (E) (see §|2J, 
we square, divide by the orbital period, and average over 
\A\ and Cl: 



<(Ai?) 2 ) 



AJ 2 (E,J,\A\,n) 



\A\,Q 



P{E)Ji{E) 



(80) 



The diffusion time scale near the loss cone boundary at 
J 2 = J 2 C (E) ~ GM, rbin, where r\>in is the radius from 
the center of the binary within which gravitational sling- 
shot can occur (see §EJ, will then be: 



((Ai?) 2 ) lc 

G 2 M 2 r 2 in P{E ) 

AJ 2 (E,J lc ,\A\,n) 



(81) 



\A\,n 



To estimate thrown we ignore the stellar self-gravity; the 
test star is then orbiting in the approximately-Keplerian 
potential of the MBH binary. Consider an elliptical orbit 
with the closest approach to the binary at r_ > rbin ■ Let 
if denote the angle between A and J (see Fig. [SJ. Since 
dt = (r 2 /J)dd where is the azimuthal angle in the 
orbital plane, A J 2 becomes: 



A T 2 

^'-'Kepler 



i-2-k 

-2|.4|sinp / r 3 sin(0)d0. 
Jo 



(82) 



The perturbation in the orbit is assumed linear in .4.; 
therefore it is sufficient to substitute in equation l|82|l 
the Keplerian orbit: 



-(0) = 



GM. 



(1-e 2 ) 



2E 1 + e cos(0 - ip) 1 



(83) 



where ip is the inclination of the orbit with respect to the 
projection of A on the orbital plane. Integration yields: 

(GM \ 3 
-j^-j eVl -e 2 stop sin V. ( 84 ) 

Note that in Kepler's potential, e 2 = 1 — R; J C (E) — 
GM./V2E; and P(E) = 2ttGM.(2E)~ 3 / 2 . Therefore, 
knowing that R\ c <C 1 at the loss cone boundary, we 
have: 



e 2 {l-e\ 



R\, 



GM. 

V 2E 



(85) 



Substituting equations H84|) and (|85|l in equation 181f) 
gives the time scale on which Brownian motion refills 
the loss cone: 



^brown, Kepler 



(GM.) 3 / 2 r hia ( GM. 
6ttA 2 { 2E 



-7/2 



(86) 
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The time scale was derived ignoring the contribution to 
the mean potential from the stars in the galaxy. If the 
galactic potential had been taken into account and the 
total potential had not been Keplerian, the functional 
dependence on E in the time scale l(86|) would have had 
a different form. The scaling in the other parameters 
(A, rbi n , M,), however, will be the same in the generic 
case. Substituting the expression for the acceleration 
from equation (|77jl yields: 



tbi. 



• 400 Myr x 



/ M. 

X \1O 6 M 



Ohard 

2-3/ 



TO* 
Wr. 



z 14/a K{x), (87) 



where ahard = GM m /8a 2 is the hard binary separation, 
a = d In M, /din a is the slope of the relation between 
the black hole mass and the velocity dispersion, and 
K[x) is a dimensionless factor encapsulating the de- 
pendence on the energy E = 2<j 2 x. Using a ~ 4.5 
ijFerrarese fc Merrittl 12000^ and rbi n /ahard ~ 0.2 in a 
galaxy like M32, we obtain ibrown ~ 5 Gyr at energies 
corresponding to the radius of influence of the black hole 
where K(x) « 1. The total mass inside the loss cone at 
these energies is small; thus the Brownian motion-driven 
diffusion of stars into the loss cone does not appear to be 
efficient in galaxies. 




Fig. 8. — Keplerian orbit of a star around the MBH binary. 
The binary is centered at the root of the angular momentum 
vector in the drawing. 

The effect of the Brownian motion in galaxies could 
be enhanced if clumpy inhomogeneities a r e pre sent in 
the galactic potential. iBacker k, SramelJ l)1999j) point 
out that the mass of inhomogeneities at the center of 
the Milky Way varies with the projected distance from 
position of Sgr A*, assuming the distance of 8.5 kpc to 
the Galactic center, as: 



m(r) ~ 1O 3 M 



pc 



1.5 



(88) 



The accelerations implied by these inhomogeneities, A ~ 
Gm(r)/r 2 « 10~ 6 cm s _1 at r = 1 pc, are comparable to 
the accelerations due to the stars that are given by equa- 
tion H77fl and amount to A* = 0.5— 3x 10~ 6 cm s _1 using 
m* = M , M. = 2.6 x 1O 6 M , and a = 75-110 km s _1 . 
The presence of compact inhomogeneities that do not 
suffer tidal truncation in the galactic potential, such as 
intermediate-mass black holes, could in principle enhance 



the amplitude of the Brownian motion beyond the pre- 
dictions of this section. 

Scaling the results of this sect ion to the simulations 
of iMilosavlievic fc Merrittl 1)2001(1 . we find that the time 
scale on which the Brownian motion fed stars into the 
loss cone amounted to ibrown ~ 1 Myr x K(x), which 
is of the order of the dynamical time. While we argued 
above (§0) that gravitational scattering into the loss cone 
was sufficient to keep it full in the iV-body simulations, 
Brownian motion probably also contributed substantially 
to the maintenance of a full loss cone in these simula- 
tions. This underscores a critical difference between the 
published TV-body simulations and real galaxies; the TV- 
body simulations are dominated by small- TV effects and 
real galaxies are not. 

8. COALESCENCE 

When the MBHs approach each other sufficiently 
closely, emission of gravity waves becomes the dominant 
mechanism extracting the binding energy. The binary 
coalesces within a time t gr when its semimajor axis a gr 
and the ecce ntricity e gr satisfy the relation derived by 
iPetersI lfl964|) : 

4 _ 64 G 3 M 1 M 2 (M 1 + M 2 )F(e gr ) t 
a gr — ^ 5 i gr , (oa; 

where F(e) has the form: 
F{e) 
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The ratio between the hard binary separation ahard = 
G/i/4cr 2 at which MBH binaries first form, and the coa- 
lescence separation a gr for circular binaries reads: 



1/4-2/a 



-1/4 



(91) 



where, again, a is the slope of the M,-a relation, and 
p = M 2 /M 1 < 1 is the ratio of the MBH masses. The 
factor of ~ 10— 100 that the MBH binary needs to decay 
on its path from initial formation to coalescence is some- 
times called the "final parsec problem." The required 
decay factor decreases with decreasing mass ratio; how- 
ever when p < 10~ 3 , the smaller MBH's host galaxy is 
likely to be severely tidally disrupted during the inspiral 
and thus deposit the black hole on an orbit larger than 
ahard- The time scale on which the smaller MBH spi- 
rals inward due to dynamical friction could the n be too 
long to p ermit the formation of a hard binary ijMerrittl 
l2000tlYull2002]l and the MBHs would not coalesce even if 
the mechanisms discussed in this paper imply that they 
should. 

We proceed to a qualitative discussion of the scenar- 
ios for the evolution of MBH binaries with masses lying 
within several ranges of values. Bright, early type galax- 
ies (My < —21) typically have shallow nuclear density 
profil es, |a'logi//dlogr| < 1, within a break radius rt, 
(e.g..lGebhard tet alJl9 961. althoug h some except ions to 
this rule have been found (e.g.. ILaine et al.l l2003) . The 
break radius is usually located well outside the dynamical 
radius of influence of the MBH. Therefore MBH binaries 
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with initial black hole masses M. > 1O 8 M are from the 
outset situated in low-density stellar cusps and the mass 
initially inside the loss cone is small compared with M m . 
In addition, the time scales for diffusion into the loss cone 
are very long (tdiff > 10 11 yr). Therefore, binaries that 
form in the mergers of galaxies with black hole masses 
above 1O 8 M have difficulty decaying much beyond the 
orbital separation of a hard binary. In this mass range, 
ihard/ftgr ~ 10, implying that these binaries would be 
long-lived. Exceptions might be merger events with ex- 
treme mass ratios, p ~ 1CV 2 , for which the required decay 
factors are smaller. If collisional loss cone refilling is very 
inefficient, the dominant loss cone replenishment mecha- 
nism might be via chaotic or centrophili c orbits if the po- 
tentia l is appreciably non-axisymmetric l)Merritt fe Poor] 
2003). Even in axisymmetric, non-spherical potentials, 
larger ste llar masses are available inside the loss cone 
l)Yull2002(l than suggested by the numerical estimates of 
the preceding sections (generalization of the results 21"|3 
to axisymmetric potentials would be straightforward but 
tedious). Furthermore, the longer a binary persists in 
an uncoalesced state, the more likely that a third MBH 
will fall in during a subseque nt merger, facilitating decay 
by three-body interactions ijMakino fc Ebisuzakil 11994: 
I Valt onenll 1 996b|) . Thus we would still expect coalescence 
to sometimes occur for high-mass binaries. 

At the other end of the mass range, M, < few x 
1O 6 M , the black holes are expected to be embedded in 
approximately isothermal stellar cusps, |dlogi//dlogr| ~ 
2. The formidable factor that the binaries need to 
cover from formation to coalescence, ahard/agr ~ 100 
is slightly larger than for more massive binaries due 
to the M.-dependence in equation (|9*T|l . But a com- 
bination of processes contribute to effective orbital de- 
cay in these galaxies. Forem ost, V-body simulations 
(Milos avlievic fc Merrittll200lD have shown that the bi- 
naries shrink by a factor of ~ 5 — 10 beyond the hard 
binary separation immediately following the formation 
of a binary. This takes place within a dynamical time 
and is a consequence of the abundance of cold (large 
binding energy) stars in the nascent loss cone to which 
the kinetic energy of the MBHs can be effectively trans- 
ferred in slingshot ejections (see jJ52J. The ejection of 
most bound stars leads to a reduction of the density pro- 
file close to the sphere of influence of the MBH. How- 
ever there are good reasons to believe that galaxies with 
M, ~ 1O 6 M are in position to repair cusp damage on a 
time scale comparable to the binary hardening time scale, 
since both time scales are proportional to the relaxation 
time. Indeed, the phase space transport due to two-body 
relaxation will ten d to rebuild a cusp p rofile p cx r~ 7 / 4 
around the MBH l|Bahcall fc Worflll976ft . Moreover, the 
nearest and the best resolved MBH nucleus — Sgr A* clus- 
ter at the center of the Milky Way — exhibits evidence for 
young, massive stars of spectral type O -B deep inside the 
radius of influen ce of the MBH (e.g.. iFiger et al.ll200nt 
Ghez et al. 2003). This suggests that even in galaxy nu- 
clei with Msg r A* ~ 4x 1O 6 M black holes, star formation 
can repopulate the central cusp in as short a time as 10 7 
yr- 

After the loss cone is empty of its original content, 
two-body relaxation lets stars diffuse into the loss region. 
Note that the binary is still a factor ~ 10 — 20 wider than 



the separation conducive to gravitational coalescence in 
a Hubble time. From equation (|ll|l . the equilibrium dif- 
fusion supplies an influx of stars at the rate: 
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and the non-equilibrium diffusion could in principle sup- 
ply stars at a rate larger by a factor of J 7 none q /J-' eq ~ 
2 — 10. Combining equations l|92|> and l|61|) . setting 
m* = M , and using A$/2cr 2 ~ 10 which is reasonable 
for isothermal galaxies with M. ~ 10~ 3 M ga i, we find: 
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where we have explictly indicated that diffusion becomes 
the dominant driver of binary hardening only after the 
immediate post-merger shrinkage to o ~ O.ldhard has 
been completed. The factor describing the enhancement 
of the loss cone flux out of equilibrium, JFnoncq/.^eq, 
is greatly uncertain and large (~ 10) when the phase- 
space is subject to episodic perturbations excited by the 
substructure inside the galaxy (orbiting giant molecular 
clouds, super-star clusters, etc.) which can restore the 
large phase-space gradients at the loss cone boundary 
(see §0J. We integrate equation (|93|l to obtain: 



0.1a ha 
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Although this result suffers from a variety of uncertain- 
ties, it suggests that for M. < few x 1O 6 M , the binary 
separation can decay by a factor of at least a few out of 
the factor of x 10 — 20 that remain toward gravitational 
coalescence. Hence we expect that stellar-dynamical pro- 
cesses might often serve to bring MBH binaries in this 
mass range to complete coalescence. 

Finally, MBHs in the intermediate range of 1O 6,5 M < 
M, < 1O 8 M live in nuclei with a range of cusp slopes, 
0.5 < \dlogv/dlogr\ < 2 lIGebhardt et alJ ll996). The 
mechanisms discussed here would facilitate coalescence 
in galaxies with steeper cusps, while in galaxies with 
shallower cusps, binaries might be shrunk to separations 
within a factor of a few of a gr , where other mechanisms 
(triaxiality, gas dissipation etc.) could plausibly bridge 
the remaining gap. 

9. SUMMARY AND DISCUSSION 

We conclude with a tabulated summary of physical 
regimes that have been discussed so far (Table [I). 

Clearly the long-term evolution of a massive black hole 
binary can be very different in different environments. 
We distinguish three characteristic regimes. 

1 . Collisional. The relaxation time is shorter than the 
lifetime of the system and the phase-space gradients at 
the edge of the loss cone are given by steady-state solu- 
tions to the Fokker-Planck equation (§EJ). The densest 
galactic nuclei may be in this regime. Resupply of the 
loss cone takes place on the time scale associated with 
scattering of stars onto eccentric orbits. Most of this 
scattering occurs in the "diffusive" (as opposed to "pin- 
hole" ) regime and the decay time of a MBH binary scales 
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Table 1. 



Form of Decay Regime 

a -1 <x t + const pinhole 

t/N a + const, (0 < a < 1) pinhole & diffusion 

t/N diffusion 

const large-./V limit 

ln(l + t/to) + const pure reejection 

VT+ Dt pinhole with reejection 



as \a/a\ ~ ~ N. In the densest galactic nuclei, colli- 
sional loss cone refilling may just be able to dr ive a MBH 
binary to coalescence in a Hubble time (e.g. [Yu| [2002). 
A^-body studies with any feasible N are guaranteed to be 
in the collisional regime, and may also be in the pinhole 
regime corresponding to an always-full loss cone, leading 
to artifically high decay rates, a -1 ~ t (§EJ. In iV-body 
studies with N < 10 6 , Brownian motion of the binary 
further enhances the decay (§01. 

2. Collisionless. The relaxation time is longer than 
the system lifetime and gravitational encounters between 
stars can be ignored. Examples are the low-density nuclei 
of bright elliptical galaxies; dark matter density cusps are 
also in this regime if the dark matter particles are much 
less massive than stars. The MBH binary quickly inter- 
acts with stars whose pericenters lie within its sphere of 
influence; in a low-density nucleus, the associated mass is 
less than that of the MBH binary and the decay tends to 
stall at a separation too large for gravity wave emission 
to be effective. However evolution can continue due to 
reejection of stars that lie within the binary's loss cone 
but have not yet escaped from the system. In the spher- 
ical geometry, reejection implies \a/a\ ~ (1 + t/to) /a, 
leading to a logarithmic dependence of binary hardness 
on time (§ EJ. Reejection in galactic nuclei may con- 
tribute a factor of ~ a few to the change in a in MBH 
binaries over a Hubble time. 

3. Intermediate. The relaxation time is of order the 
age of the system. While gravitational encounters con- 
tribute to the re-population of the loss cone, not enough 
time has elapsed for the phase space distribution to have 
reached a collisional steady state. We argued that most 
galactic nuclei are likely to be in this regime. The flux 
of stars into the loss cone can be substantially higher 
than predicted by the steady-state theory, due to strong 
gradients in the phase space density near the loss cone 
boundary produced when the MBH initially formed (§BJ. 
This transitory enhancement would be most important 
in a nucleus that continues to experience mergers or in- 
fall, in such a way that the loss cone repeatedly returns 
to an unrelaxed state with its associated steep gradients. 

The evolution of a real MBH binary may reflect a com- 
bination of the different regimes summarized above, as 
well as other mechanisms that we have discussed only 
briefly or not at all. We note a close parallel between the 
"final parsec problem" and the problem of quasar fuel- 
ing: both requre that of order 10 8 M Q be supplied to the 
inner parsec of a galaxy in a time shorter than the age of 
the universe. Nature clearly accomplishes this in the case 
of quasars, probab ly through gas flows driven by torque s 
from stellar bars (Shlosman, Begelman & Frank 1990). 
The same inflow of gas could contribute to the decay of 



a MBH binary in a number of ways: by leading to the 
renewed formation of stars which subsequently interact 
with the binary; by inducing torq ues which extract angu- 
lar m omentum from the binary ijArmitage fe Nataraianl 
2002); through accretion, increasing the masses of one or 
both of the MBHs and reducing their separation; etc. 

In massive elliptical galaxies or spheroids, the frac- 
tion of mass in the form of gas during the last 
major merger event is like ly to have been small 
ijKauffmann fc Hae hnelt 200j3). In addition to the mech- 
anisms discussed in this paper, decay of a MBH bi- 
nary could be enhanced by the presence of a signif- 
icant population of "ccntrophilic" orbits, orbits that 
pass near the center of the potential once per crossing 
time. Such orbits may constit ute a large fraction o f 
the mass of a triaxial nucleus ijPoon fc Merrittl 12002'). 
Even in galaxies where none of these mechanisms is ef- 
fective and the decay stalls, infall of a third MBH fol- 
lowing a merger or accretion eve nt could accelerate the 
decay via the Kozai mechanism l)Blaes. Lee fc Socrates! 
120021) . or the strong three-body gravitational interactions 
ijMakino fc Ebisuzakilfl994t lValtonenllT996b|) . which in- 
duce large changes in the binary's eccentricity and en- 
hanced rates of gravity wave emission at pericenter. 

We have not attempted in this paper to make quan- 
titative estimates of the decay rates or stalling radii of 
MBH binaries in real galaxies. Such a program would be 
problematic for a number of reasons, most importantly 
the uncertain history of the loss cone and the unknown 
elapsed time and character of the most rece nt merger 
(§||}. However we argued that a recent study l)Yull2002]) 
could have overestimated stalling radii, due both to the 
neglect of the various new mechanisms discussed here, 
and also by ignoring the influence that MBH formation 
must have had on the form of the nuclear density profile: 
this profile must have been steeper before the MBH bi- 
nary formed, leading to substantially enhanced loss cone 
fluxes and more rapid decay early in the life of the bi- 
nary. Nevertheless, none of our results would lead us to 
rule out the existence of uncoalesced MBH binaries in 
at least some galaxies, particularly galaxies with large, 
low-density nuclei and a little gas (§0- 

A striking and disappointing conclusion of this study 
is the difficulty of using numerical iV-body simulations 
to follow the long-term evolution of MBH binaries (§EJ. 
We argued that a variety of collisional effects associated 
with small N would lead to a more rapid decay of the 
binary than in real galaxies. While there is still a useful 
role for iV-body simulations in following the initial for- 
mation and early evolution of a MBH binary, when much 
of the decay and associated cusp destruction takes place, 
the rapid decay of a(t) seen in most published simula- 
tions is due to loss-cone repopulation occuring at rates 
much higher than expected in real galaxies. We believe 
that the future role of ./V-body simulations in this field 
will be limited to predicting the form of the phase-space 
gradients produced during the formation of a MBH bi- 
nary, which can then be adapted as initial conditions for 
a collisionless or Fokker-Planck treatment. A carefully 
designed set of iV-body experiments might also be used 
to understand the TV-dependence of collisional loss cone 
effects which could then be extrapolated to the larger- N 
regime characteristic of real galaxies. 
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